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An efficient tocreaenteHterntive approach Sir differentiating ad vaaced Row cocks is successful? doaonstrated 
on a two-dimensional kvlsdd Dodd problem. The merited employs the reverse-mode capability of the automatic- 
diflerentiation software tool ADIFOR 3.0 and is proves to yield accurate first-order aerodynamic sensitivity 
derivatives. A substantial reduction in CPU time and computer memory Is demonstrated to comparison with 
results from a straightforward, black-box reverse-mode application of ADIFOR 3.0 to the same flow code. An 
ADIFOR-assfated procedure for accurate second-order aerodynamic sensitivity derivatives is successfully verified 
on an tavfadd tr an soni c lifting airfoil example problem. The method requires that first-order derivatives are 
calculated first using both the forward (direct) and reverse (adjoint) procedures; then, a very efficient noniterative 
calculation of afi second-order derivatives can be accomplished. Accurate second derivatives (Le^ the complete 
Hessian matrices) of lift, wan drag, and pitching-moment coefficients are calculated with respect to geometric 
shape, angle of attack, and freestream Mach number. 


L Introduction 

T HIS paper revisits and focuses entirely on the computational 
challenges that are associated with the efficient calculation of 
aerodynamic sensitivity derivatives (SDs) from advanced computa- 
tional fluid dynamics (CFD) codes. Of course, an accurate efficient 
methodology for obtaining these derivatives is a critical prerequi- 
site concern that must be addressed first by the aerodynamic design 
engineer who chooses any gradient- based methodfs) for design op- 
timization and/or for estimating quantities related to aerodynamic 
uncertainty. Thus computing SDs from high-fidelity nonlinear CFD 
codes is an enabling technology for design of advanced concept 
vehicles. 

In recent years significant progress has been achieved in the ef- 
ficient calculation of accurate SDs from these CFD codes. 1 The 
automatic differentiation (AD) software tool ADIFOR (Automatic 
Differentiation of FORTRAN) has beat proven an effective tool far 
extracting aerodynamic SDs from these modem CFD codes. 2-6 The 
foundation of the present work is found in Refs. 3 and 6; the present 
study builds on these earlier studies in an effort to exploit the full 
potential of the latest version of ADIFOR 3.0 (Ref. 7) for obtaining 
SDs from CFD codes. 

In Ref. 2 a strategy known as the ADD method was first pro- 
posed and lata- successfully demonstrated in Ref. 3, whereby AD 
was applied to a CFD code in incremental-iterative (I-I) form. The 
ADI1 method is a hybrid (compromise) scheme, designed to main- 
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tain as much as possible the computational efficiency of a hand- 
differentiated (HD) approach and the ease of implementation of 
a straightforward black-box (BB) application of AD; at the same 
time the accuracy of the SDs is not compromised. A comprehensive 
overview of the development of the ADD scheme is given in Ref. 
3. Also included in Ref. 3 is a comparison of die ADII scheme with 
the HD and BB approaches; computational issues associated with 
CPU time, computer memory, and SD accuracy are discussed. The 
two-dimensional effort of Ref. 3 was later extended to the three- 
dimensional code CFL3D, including “in-pantileT computation of 
the derivatives.®- 9 Appropriate references to die version of CFL3D 
used can be found in Ref. 6. 

The success reported in these previous works 3 - 8,9 could be con- 
sidered limited, however, because all ADIFOR implementations 
reported therein were forward-mode (direct) differe ntiati ons. It is 
very difficult to make any forward-mode implementation of deriva- 
tive calculations computationally competitive with a reverse-mode 
(adjoint) implementation whenever the number of design variables 
(NDV) of interest is considerably larger than the number of out- 
put functions (NOF) of interest, and NDV much greater than NOF 
is more typical for aerodynamic design problems. In recent stud- 
ies the new reverse-mode capability of ADIFOR 3.0 (not available 
for the earlier referenced studies) has been successfully verified in 
Ref. 6 by application to a parallel version of CFL3D and in Ref. 
10 by application to a sequential linear aerodynamics code. These 
applications resulted in accurate design SDs as well as stability and 
control derivatives, respectively. The application reported in Ref. 6 
involved BB AD of the entire CFD code, but iterative execution of 
the reverse mode was required only over die last iteration of the 
function evaluation. 

In the present study it is proposed and demonstrated that the 
reverse-mode capability of ADIFOR 3.0 can also be applied to CFD 
codes in I-I form, resulting in a hybrid adjoint-variable (AV) scheme 
(known herein as the ADII-AV method) that is analogous to the 
forward-mode ADII scheme of Ref. 3 and elsewhere. The motiva- 
tion of this new reverse-mode ADII-AV scheme is identical to that 
of the earlier forward-mode ADII method: greater computational 
efficiency is sought over a BB implementation of AD, without any 
loss of accuracy in the calculated SDs and without unmanageable 
complications upon implementation. 
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Following development of the proposed new ADII-AV scheme, 
the second focus of the present study is that of calculating second- 
order (SO) aerodynamic SDs from CFD codes. The motivation for 
calculating SO SDs is to advance the possibility of (or greater capa- 
bility of) 1) second -order gradient-based aerodynamic design op- 
timization, 2) analysis and design involving vehicle stability and 
control and 3) robust design [i.e., design under uncertainty, where 
a first-order, second-moment method requires SO SDs 11 ]. This sec- 
ond part of the present study is another extension of Ref. 3 , wherein 
the computational issues associated with calculating these higher- 
order derivatives were addressed, and sample calculations of SO 
derivatives using AD were reported from a two-dimensional CFD 
code. 

In Ref. 3, four procedures for calculating SO CFD SDs were pro- 
posed, but only one of the less efficient methods was actually tested; 
ADIFOR 3.0 currently provides three forward-mode variations Ra- 
the cdculation of SO SDs by similariy inefficient methods. The most 
efficient (for large NDV) SO SD scheme was not tested in the earlier 
study 3 but has been successfully implemented in the present study. 
Reverse-mode (adjoint-based) differentiation is required within this 
efficient SO SD scheme, via either HD or AD. However, with the 
availability of ADIFOR 3.0 and the new ADII-AV scheme the door 
has been opened for AD implementation and testing of this SO SD 
scheme for CFD codes. The results of this effort to date are reported 
here. These efficiently computed SO SDs have been used to demon- 
strate an approach for CFD input uncertainty propagation and robust 
design optimization fen- a quasi-cme-dimensional flow application in 
Refill. 

II. Basic Equations and Theoretical Development 

The equations summarized subsequently are discussed in greater 
detail in the references, in particular. Ref. 3. These concepts are 
known in the mathematical optimization community, 12 but the de- 
tails developed here do not appear to be generally known throughout 
the CFD community. The aerodynamic output functions of interest 
F and the discretized conservation laws of steady compressible fluid 
flow R, including boundary conditions, can be represented symbol- 
ically as follows. 

Aerodynamic output functions: 


In the preceding equations i, j, and / are “free” indices, and re- 
peated indices m and p are (by convention) “summation” indices. 
The reverse-mode (adjoint) approach for the FO SDs is developed 
starting with an application of the chain rule: 


3Fj 3 R, _ 3Fi_ 
3Ri BQm ~ BQm 


in order to show that the conventional adjoint variable is indeed a 
derivative. With Eq. (5) it then follows from Eq. (4) that Eq. (3) can 
be written as 


f'= a Il 

" dbj 


8R, pi dbj) 3X„ ri dbj 


( 6 ) 


A more conventional derivation of the adjoint-variable method (i.e., 
the Lagrange-multiplier method) gives 


F' = 
r U 


dbj 


j (3R, r 


SKA 

dbj) 


+Hl x > 

+ ax P PJ+ dbj 


(7) 


3 R, a Ft _ 0 

Gim = h, W„ + W m -‘ m 


( 8 ) 


where A.,-/ is called the adjoint variable. Comparison of Eqs. (6) and 
(5) with Eqs. (7) and (8), respectively, reveals the identity 


3Ri 


— Ft 


(9) 


One objective of this particular development of the AV method for 
aerodynamic SDs is to ensure that the relationship given by Eq. (9) 
is clearly understood; that is, A.« is the derivative of the output F, 
with respect to the intermediate variable Ri and is accumulated in 
the reverse-mode AD. 

The I-I strategies for solving the preceding equations for Q' mj 
and/or X it required for F' tj are reviewed here; additional detail is 
found in Refs. 3 and 13 and elsewhere. The I-I method for solving 
the nonlinear flow of Eq. (2) is 


F = F\Q(b),X(b),b] (1) 


2 N+ I _ n N pN pN 
m — Um ~ “ ml K t 


( 10 ) 


Nonlinear state equations: 

* = R[0(»),.XW,&] = O (2) 

where Q is the vector of state (field) variables, X is the vector of 
computational grid coordinates, and b is the vector of input (design) 
variables. 

A. First-Order Sensitivity Derivatives 

Subject to the following definitions, index (summation) notation 
is now introduced: 



(This notation will be necessary to avoid subsequent ambiguity when 
the SO SD methods are presented.) The forward-mode (direct) ap- 
proach for calculating first-order (FO) SDs is developed by differ- 
entiation of Eqs. (1) and (2) with respect to the design variables; the 
result is 


d F t 3F, 3 Fj 3Fi 

>J dbj ~ 3 Q m VmJ 3X P PJ dbj 


_ d R, _ 3R, 
,J ~ dbj~ 3 Q m 


a r, 
ax. 




3 R, 
3bj 


(3) 

(4) 


where the superscript N is the iteration (pseudo-time step) index 
and the operator 


pN _ 

* 3 R? 



(ID 


represents the solution algorithm of the particular CFD code of 
choice. The tilde in Eq. (11) serves to indicate that P% can be viewed 
as any computationally efficient approximation (often a very crude 
approximation) of the exact operator associated with true Newton- 
Raphson iteration. Thus the CFD solution algorithm is simply quasi- 
Newton iteration. 

The I-I method for solving the forward mode FO SD equation (4) 
is 


O'" 

V£ mj 


H + l _ n M p p,U 
— Qmi~ "mlR ll 


where superscript M is the FO SD iteration index and 


R 


,m _ %Ri 
“ ~ BQm 


■ mj 


+ 


SRlyl , &Rl 

dT p pj+ 3b] 


( 12 ) 


(13) 


With the I-I methodology the CFD flow solution operator P„i is 
also used to solve the SD equations; this operator in Eq. (12) is 
evaluated and fixed using the steady-state solution for the nonlinear 
flow. The requisite terms of Eq. (13) are constructed either by hand 
differentiation (i.e., the HDII method, which is very tedious and time 
consuming to complete with accuracy for advanced CFD codes) or 
by AD, which is the forward-mode ADD method of previous studies. 
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In contrast with the ADB method, a strai gh t fo rw ar d BB appli- 
cation of AD to the CFD code, which is the ADBB method, is 
represented symbolically as 


w = &Ft r , _dGr»_/dG\ 

db k Wm 3M«/’ " ibt \d b)^ 


C = <-*rC*f (!4) 

Oeatiy ADB [Eqs. (12) and (13)] and ADBB [Eq. (14)] yield the 
same result at steady-state convergence of each [recall Eq. (2)]; 
however, ADB is potentially more efficient than ADBB because of 
user intervention in the application of AD. With ADB the following 
is trim: 

1) The operator P^ can be evaluated only once [hence denoted 
Pmi in Eq. (12)] using the steady-state field variables Q and then 
raised for all M iterations and for all j = NDV design variables in 
obtaining the Q' nj . 

2) All derivatives except Q' mJ cm be computed once outside the 
iteration loop and frozen for reuse inside die loop. 

3) Evaluation of the terms P’^j in Eq. (14) can be avoided com- 
pletely for all iterations and all design variables. 

The I-I method for solving the reverse-mode AV, FO SD Eq. (9) 
is 


where 


i*+i _ i* 
a i 7 ~ K il 


P~iG, 



3F, 

SQm 


(15) 

(16) 


The requisite terms of Eq. ( 1 6) are constructed either by hand (Le., 
the HDII-AV method, having the same drawbacks as the forward- 
mode HDB method) or by AD, which is the proposed new ADB- AV 
scheme. The BB AD in reverse mode (the ADBB-AV method) has 
been verified in Ref. 6. The objective of the proposed ADI1-AV 
scheme is improved computational efficiency over the ADBB-AV 
approach without resulting loss of accuracy or significant loss in 
the ease of implementation. The mechanisms from which improved 
computational efficiency can be expected are analogous to those 
explained before when the forward-mode ADD and ADBB methods 
were contrasted F ur t h ermore, the ADB-AV scheme should lend 
itself to more permanent generalized coding implementations than 
the ADBB-AV approach. This is because with the ADB-AV method 
the manna' in which AD is applied is independent of, yet valid for, 
all of the particular aerodynamic inputs and outputs of interest 
The forward-mode application of ADIFOR produces FORTRAN 
source code for very efficient calculation of the vector (or ma- 
trix) product that results horn the postmultiplication of a large 
Jacobian matrix by a known input vector (or matrix). This at- 
tribute of forward-mode AD is exactly what was requited to con- 
struct the ADB method; specifically, file terms (dfii/dQ m )Q'* and 
(dRi/3X p )X' P j of Eq. (13) are of this type. In contrast, however, 
the forward-mode application of ADIFOR produces source code 
that is prohibitively inefficient for calculating tbe premultiplication 
of a large Jacobian matrix by a known input vector (or matrix). This 
weakness of the forward-mode application of ADIFOR is exaedy 
the strength of the reverse-mode option now available in ADIFOR 
3.0. Thus, the proposed new efficient ADB-AV scheme has become 
possible with this reverse-mode capability. That is, through reverse- 
mode application of ADIFOR 3.0 it is now feasible to produce (au- 
tomatically) the source code required for efficient evaluation of the 
tenn Xjf ( 3Ri/dQ m ) in Eq. (16), that is, the premultiplication of a 
large Jacobian matrix by a vector. 


B. Second-Order Sensitivity Derivatives 
The SO SD methods are presented in index notation subject to 


the following definitions: 


_ #F, /d>F\ 

ijk ~db k dbj VdW,,/ 

_ d 2 r, (#r\ 

,ik db k dbj \<&) lJt 

„ _ cPe* _(#q\ 

mJk *t> k dbj \dip) Kjk 

_ d*X p (<fx\ 

Pjk db k dbj V<*V W * 


The following differential operator is also introduced for subsequent 
notational compactness: 


80-, +30.W , £0 
D b* ag« 3X» 3b* 


( 17 ) 


where repeated indices n and q are summation indices. 

Differentiation of the FO forward-mode equations (3) and (4) 
with respect to the design variables yields SO method 1: 


F „ = 

iJt db*db y 
d 2 R, a R, 


ELq» +^LX" 

3Q m Umjk ax p pjk D b k 


(18) 


R " ik ~ 6b, dbj 




o" I aRl X" 4- 
T- A *-t- 


ax 


dr;. 


= 0 ijt (19) 


The terms of DF,'/Db* and DR^/Db* are many and very com- 
plicated; detailed expansion of these toms is provided in tire Ap- 
pendix. Using symmetry of the Hessian Q" yt = Q'^j, SO method 1 
requires (NDV 2 + NDV)/2 solutions of the large linear systems 
of Eq. (19) for 2" yt ; in addition, the method requires NDV solu- 
tions of Eq. (4) for the FO SDs Q' mj . SO method 1 was verified for 
a two-dimensional CFD code in Ref. 3 by ADBB differentiation of 
tbe code’s existing HDB scheme [Eqs. (12) and (13)] for tbe FO 
SDs. 

Alternatively, differentiation of tbe FO reverse-mode, Eqs. (7-9), 
with respect to the design variables, yields SO method 2: 


F" — 
t ijk ~ 


d 2 F, „ /3R., ^3RA 

dbtdbj l!l \ax p PJ+ dbj) 


[a Fi 

aR,\ 

i „ DF u 


+ W 

+Xil aT p ) 


(20) 

dG im 

„ 3R/ 

DG,_ 

(21) 

imk ~ db* 

~ i,l aQ m 

+ 

2 
^ i 

ii 

Q 

£ 


SO method 2 requires NDV x NOF solutions of the large linear sys- 
tems of Eq. (21) for X' ilk ; in addition, the method requires NDV 
solutions of the FO equation (4) for Q' mJ plus NOF solutions of the 
FO equation (9) for A.,/ . This SO method 2 is eliminated from further 
consideration because it is unconditionally less computationally ef- 
ficient than the remaining two SO SD methods. 

Introduction of tbe AV approach within SO method 1 to eliminate 
Q'mjk yields SO method 3: 




( m 




+ 


Efw., 

Db* + “ Db k 


( 22 ) 


SO method 4 is similar and computationally equivalent to SO 
method 3 and is developed by introduction of the AV approach 
within SO method 2 to eliminate k' llt ; the result is the identity 


^il 


£ 

Db* 


Q 


, D Gun 
mj Dbt 


(23) 


where SO method 4 uses Eq. (23) to replace equivalent toms within 
SO method 3. The equivalent SO SD methods 3 and 4 do not re- 
quire solution of large systems of linear equations for higher-order 
derivatives such as Q" mjl or k' ilk . These two SO SD schemes do, 
however, require solution of both forward-mode and reverse-mode 
equations (4) and (8) for Q' mj and A. ( j, respectively. This is a total of 
only NDV +NOF solutions of large systems of linear equations. 

One significant conclusion of the preceding analysis is that SO 
methods 3 or 4 should be computationally more efficient when- 
ever NDV 2 +NDV is greater than 2 x NOF. With typical design 
problems in aerodynamics, NDV is often much larger than NOF; 
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typically, NOF is three or less (often only one), whereas NDV is on 
the order of tens to hundreds. The advantage in favor of method 3 
ot 4 for SO SDs is then overwhelming because of the NDV 2 term, 
which dominates. Once both the forward-mode and reverse-mode 
schemes are in place for calculating the FO SDs, then complete 
SO SD information is available almost “for free,” that is, the SO 
SD are obtained through an explicit, noniterative calculation. The 
source code for implementation of method 3 or 4 is constructed 
“automatically” via BB application of the forward-mode capability 
of ADIFOR to appropriate pieces of the existing source code from 
which the FO SDs are obtained. For example, the extremely com- 
plex terms DFf/Db*, DF^/Db*, and/or Dtj lm /Db t (see Appendix) 
of methods 3 and/or 4 are easily constructed with a forward-mode 
application of AD. 

In Ref. 3, SO methods 3 and 4 were proposed but not actually 
tested. Consequently, oik primary goal erf the present study is suc- 
cessful implementation and verification of the highly efficient SO 
method 3 (or equivalently, method 4); method 3 is actually chosen 
in this study. 

ILL Results and Discussion 

A. First-Order Sensitivity D eri v a ti ves, 

ADII-AV Method, Model Problem 

The proposed ADII-AV method [Eqs. (15) and (16)) has been 
successfully implemented in a CFD code and verified for accuracy 
on a simple two-dimensional inviscid internal flow model problem. 
This CFD code solves the two-dimensional Euler equations by a 
conventional upwind finite volume approach on a very coarse grid 
but one that is sufficient for computationally verifying SDs. As ex- 
pected, when FO SDs computed by the new ADII-AV scheme are 
compared with SDs computed by a hand-differentiated implemen- 
tation of Eqs. (15) and (16) (Le., the HDII-AV approach), the results 
are the same at convergence, as well as at each 1-1 step. In addition 
the accuracy of the computed SDs has been successfully verified by 
a finite difference method. 

Preliminary tunings were conducted on a Sun workstation to eval- 
uate the potential for improved computational efficiency of the new 
ADH-AV scheme with respect to the ADBB-AV approach of Ref. 6. 
Computational tuning comparisons are given in Thble 1, which fo- 
cuses exclusively on AD performance. Therefore, relative tunings 
are given as CPU time per iteration per grid point per differentiated- 
aerodynamic -output function. Furthermore, each timing result has 
been scaled by the comparable timing result obtained from the very 
efficient hand-differentiated reverse-mode scheme (i.e., the HDII- 
AV method). Table 1 illustrates that although the new ADH-AV 
scheme is almost five times slower than the efficient HDII-AV 
scheme it represents a substantial improvement over results obtained 
from the straightforward black-box procedure (i.e., ADBB-AV is 
about eight times slower than HDH-AV). 

The improvement in computational efficiency achieved to date 
is substantial when the reverse-mode application of ADIFOR 3.0 
in incremental-iterative form is compared with the black-box ap- 
proach. Furthermore, the timing result for the ADH-AV scheme is 
projected to improve by an additional 30% over that repotted here. 
Thus the relative timing given in Table 1 for ADII-AV /HDII-AV is 
projected to drop from 4.7 to about 3.3. This projection is based on 
using a strategy where the forward-pass execution of the ADIFOR- 
enhanced, reverse-mode code will be performed only once (instead 
of during each iteration) in order to create the required ADIFOR log 

Table 1 Relative CPU timing 
comparison, model problem 


files. Thereafter, by repeatedly reusing these fixed log files only re- 
verse passes will be repeatedly executed during the iterative solution 
process for all aerodynamic output functions of interest. 

Another important computational concern mitigated by the new 
ADH-AV method is computer memory, particularly the issue of large 
disk files created during execution of reverse-mode derivative code 
created by ADIFOR 3.0. With the black-box (ADBB-AV) approach 
these large ADIFOR log files (which are created on a forward-pass 
execution and are read during die reverse pass) will accumulate and 
become larger with every iteration of the ADIFOR-enhanced flow 
code. This file growth can rapidly deplete the available disk space, 
even on the largest computers. In Ref. 6 this difficulty was addressed 
by development of the iterated reverse-mode scheme, where only 
the log files for the final forward-pass iteration are stored and used 
during the subsequent iterative solution for the derivatives. With 
the ADH-AV approach, however, the required disk space is not as 
restrictive an issue because it remains fixed and does not accumulate 
during the iterative solution process. In the present example the total 
storage requirement for log files with the ADH-AV method is only 
64% of that required for a single iteration of the ADBB-AV method. 

In addition to the log-file disk memory, required only for the new 
reverse-mode capability of ADIFOR, there are substantial additional 
core memory requirements. For the forward mode the core memory 
increase of the AD-enhanced code is approximately NDV times the 
core memory requirement of the original (undifferentiated) code. 
For the reverse mode the corresponding increase is NOF time that 
of the original code. 

B. Second-Order Sensitivity Derivatives, 

SO Method 3, Airfoil Example 

Results are presented subsequently from the successful verifica- 
tion of the proposed efficient noniterative SO method 3 [Eq. (22)] 
for computing SO SDs. The example problem is steady transonic in- 
viscid flow over a NACA 001 2 airfoil with freestream Mach number 
Afoo 0.80 and angle of attack a 1 deg. The two-dimensional Euler 
equations are solved on a Sun workstation in double precision using 
a conventional finite volume upwind flux- vector-splitting scheme. A 
C-mesh computational grid is used with dimensions 129 x 33 grid 
points. High-quality lift-corrected boundary conditions are used at 
the far fiefd boundary, which is placed approximately five chord 
lengths from the surface of the airfoil. 

In the present example derivatives of three aerodynamic out- 
put functions are considered; CL, C d , and Cm (Le., coefficients of 
lift, wave drag, and pitching moment, respectively). The computed 
steady-state values of these aerodynamic force coefficients are given 
in Ikble 2. Note that the number of digits shown in Table 2 (and also 
Tables 3 and 4) is to illustrate consistency rather than accuracy. In 
addition, derivatives with respect to three aerodynamic input vari- 
ables are considered. They are g (a geometric-shape variable), a, and 
Afoo- The geometric-shape variable g is a single arbitrarily selected 
y coordinate of the computational grid on the surface of the airfoil: 
simple, but one that is sufficient for verifying geometric- shape SDs. 

Calculation erf SO SDs by SO method 3 requires that all FO 
SDs are calculated first using both the forward-mode [Eqs. (3) and 
(4)] and the reverse-mode [Eqs. (7) and (8)] approaches. The cal- 
culated FO SDs from a hand-differentiated incremental-iterative 
(HDH) implementation of these two approaches are presented in 
Table 3, where the results are seen to agree, as expected. The FO 
SDs presented in Table 3 have been thoroughly verified for consis- 
tency through a meticulous implementation of the method of central 
finite differences, where agreement to six significant digits or greater 
is noted in all comparisons. 


Reverse-mode 
method tested 

Relative 

tuning* 

Table 2 Aerodynamic force 
coefficients, airfoil example 

HDD-AV 

ADII-AV 

1.0 

4.7 

Coefficient 

Value 

ADBB-AV 

7.9 

Cl 

Cd 

+02830659E+00 

+02070493F-01 

CPU tmWiaHitkjQ/ghd^atf/oitfput 
function. 

Cm 

— 02876639E— 01 



1228 


TAYLOR ET AL 


Table 3 First-order sensitivity derivatives, 
HDn method, airfoUesMapjc 


Derivative bj 

del I~ 

*>j “ 

Mqo 

dC p x 

dfi, “ 

Woo 

dfi; “ 


Forward mode 
+0.14Q5406F+00 
— 0.1087323£— 01 
— 0.4672729£ —0 1 
+0.1761807£+02 
+0.1 158625E+01 
— 0.2382580E+01 
+O.3171492£+Ol 
+03955598E+01 
— 0.16030Q2£+01 


Reverse mode 

+0. 1405406E +00 
—0. 1087323£— 01 
-0.4672729£ —01 
+0.1761807£+02 
+0.1158625£+01 
— 0.2382580£ +01 
+0.3171492£+01 
+03955598£ +01 
— 0.1603002£+01 


Table 4 Secoad-order s twri tiv tty derivatives, 
SO method 3, airfoil example 


Derivative bj/b k g 


Mo, 


<?C L g 

Ujdbk J 

Moo 

(PCp g 

dbjdbt “ 

Mqo 

&Cu g 

ibjibk “ 


+O.248807£+O3 
+0.25G402E +03 
+0.2Q5825E +03 
+0.71777E+02 
+0.134379E+02 
+0.101304E+02 
— 0.663590£+02 
-0.6024415+02 
— 0.490399£ +02 


-+0.250402E+03 +0.205825E+03 
+O.184277E+05 +0.160858E+05 
+0.160858E+05 +O.133087E+O5 
+0.134379E+02 +0.101304E+02 
+0.959310£-+03 +0.804021 £+03 
+0.804021E+03 +0.662088E+03 
-0.602441 £+02 -0.490399E+02 
— 0.449 198£ +04 -0386943+04 
— 0386943£+04 -0.3205 12£ +04 


Tables Relative CPU timings: complete SO 
method 3, airfoil example 


Computational procedure 

% of total 

Nonlinear flow, Eqs. (1) and (2) 

5.4 

Forward-mode FO SDs, Eqs. (3) and (4) 

253 

Reverse-mode FO SDs, Eqs. (7) and (9) 

69.0 

SO SDs, Eq. (22) 

0.1 

Total 

100.0 


Table 5 shows that using an ADIFOR-assisted second differentia- 
tion SO SDs can be obtained extremely fast, if one already has both 
the forward-mode and reverse-mode FO SDs. 


IV. Conclusions 

An efficient incremental-iterative approach for differentiating ad- 
vanced CFD flow codes has been successfully demonstrated on a 
two-dimensional inviscid model problem. The method employs the 
reverse-mode capability of the automatk-differentiation software 
tool ADIFOR 3.0 and has been shown to yield consistent first-order 
aerodynamic sensitivity derivatives. A substantial reduction in CPU 
time and computer memory has been demonstrated by comparison 
with results from a straightforward, black-box reverse-mode appli- 
cation of ADIFOR 3.0 to the same flow code. 

A computationally efficient ADIFOR-assisted procedure for con- 
sistent second-order aerodynamic sensitivity derivatives has been 
successfully verified on an inviscid transonic lifting airfoil exam- 
ple problem. Accurate second derivatives (i.e,, the complete Hes- 
sian matrices) of lift, wave drag, and pitching-moment coefficients 
with respect to geometric shape, angle of attack, and freestrcam 
Mach number have been calculated. With the present procedure 
second-order derivatives are now computationally feasible, at least 
in two dimensions. The computation of second+xder derivatives 
in three dimensions appears to be within reach, bin remains to be 
investigated. 

This second-order method requires that first-order derivatives be 
calculated using both the forward (direct) and reverse (adjoint) pro- 
cedures; then, second-order derivatives can be obtained in a noniter- 
ative calculation that is computationally very efficient An ADIFOR 
differentiation is used to generate a number of required second -crxler 
terms in this noniterative calculation. If one already has either to- 
ward (NDV solutions) or reverse (NOF solutions) FO SDs, then 
upon obtaining the other FO SDs (NOF or NDV additional solu- 
tions, respectively) one calculates all of the SO SDs (NOF x NDV 2 
derivatives) very efficiently. 


The SO method 3 is implemented by application (in the to- 
ward mode) of ADIFOR to appropriate pieces of the FORTRAN 
code used earlier for hand-differentiated forward-mode calculation 
of the FO SDs. The calculated SO SDs from this implementation 
of SO method 3 are presented in Table 4. The SO SDs of Table 
4 have been thoroughly verified for consistency through a meticu- 
lous application of central finite differences applied to FO SDs ob- 
tained by the hand-differentiated methods already described. Agree- 
ment to five significant digits or better is noted in this verification 
study for all SO SDs reported in Thble 4. This verification study 
was not conducted using finite differences applied to the origi- 
nal nonlinear flow code; that approach has been documented to 
be vulnerable to severe numerical inaccuracy when SO SDs are 
calculated. 3 The symmetry of the calculated SO SDs shown in Ta- 
ble 4 is expected and results from the computations performed, 
that is, no derivative symmetry was explicitly imposed on the 
problem. 

For the present airfoil example problem Table 5 illustrates (in 
terms of percentages of the total) the breakdown of relative CPU 
timings to the computational steps of SO method 3 procedure for 
calculating the SO SDs. Not included in Table 5 is the CPU time 
to the grid generation and the grid-sensitivity derivatives, negligi- 
ble for this particular two-dimensional example. Thbie 5 illustrates 
clearly the computational efficiency of the SO method 3 to SO SDs. 
Recall that results of the present example are for three aerodynamic 
output functions and three input (design) variables, where the com- 
putational work of the forward-mode and reverse-mode procedures 
for FO SDs should be approximately equal (in theory, to hand- 
differentiated code, as used hoe). In this example, however, Thbie 5 
reveals that the reverse mode was much more costly than the for- 
ward mode; apparently the three linear systems for the reverse mode 
are stiffer than the three to the forward mode. This characteristic 
of the adjoint equations has been observed by others. 7 As expected. 


Appendix: Expansion of Terms 

In this Appendix the terms DF,' ; /Di>t, DRJ ./Dfit, and DG, m /Dfit 
are expanded using the index notation already established. The ex- 
pansion of DFJj/Dbt is 


aF !j 


E = _ 

D h t dQ n BX q ~ r db k 

d 2 Ft 


96.3(2, 

d 2 Fi 


a 2p. a 2 F, 

Q'A + wir x '* Q « + 


dQ n usi p 

d 2 Fi 


dQ n dbj 

a 2 Fi 


+ ax,ae„ QiliXqk + ax q dx Xp>Xqk + ax q ab Xqk 

9 2 F. . B 2 Fi , d 2 F t 


+ "■ ~ ~ Qm ‘ + ab k ax Xpl + db k abj 


db k a<2. 


(Al) 


In Eq. (Al), the indices i, j, and k are free and repeated indices n, 
m, p, and q are summed. The toms of D/fy /Dfi* are obtained from 
Eq. (Al) by replacing everywhere F{j with R' tj and F,- with Ri (and 
thus / replaces i as a free index in the resulting expressions). 
Finally, the expansion of the toms for DG im /Db k is 


DGiw 3 Gim , dG im dG lm 
Db t ~ 3fi, Umt + 3X q qk d b k 


— hi 


a 2 r, 
3Q.3 Q, 


Q'ok + 


* Ft „ , , , a 2 */ r 


3 2 F, a 2 R, a 2 F t 

BX q dQ m " 31*36,, db k dQ m 


(A2) 
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In Eq. (A2) the indices i, m, and k are free, and repeated indices l, 
n, and q are summed. 
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